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Supermassive black holes (SMBHs; M > 10 5 M ) are known to exist at the 
centre of most galaxies with sufficient stellar mass. In the local Universe, it is 
possible to infer their properties from the surrounding stars or gas. However, 
at high redshifts we require active, continuous accretion to infer the presence 
of the SMBHs, often coming in the form of long-term accretion in active galac- 
tic nuclei. SMBHs can also capture and tidally disrupt stars orbiting nearby, 
resulting in bright flares from otherwise quiescent black holes. Here, we re- 
port on a ~ 200-s X-ray quasi-periodicity around a previously dormant SMBH 
located in the centre of a galaxy at redshift z = 0.3534. This result may open 
the possibility of probing general relativity beyond our local Universe. 

Tidal disruption of stars as a means to fuel active SMBHs was originally proposed in 
1975 (1) but it was not until over a decade later that the possibility of using the expected elec- 
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tromagnetic flares to study non-active SMBHs was proposed (2). Since then, various candidate 
tidal disruption flares (TDFs) have been identified based on luminous flares observed from opti- 
cal to X-rays (3-7), but it is only in recent years that such objects have been confirmed through 
the observation of relativistic flares (8-12). 

Swift J 164449.3+573451 (hereafter Sw J1644+57) was detected (8) by the Swift Burst Alert 
Telescope (BAT, 15-150keV) on 28 March 2011 as it reached X-ray luminosities greater than 
~ 10 48 erg s _1 . Prompt, multi-wavelength observations spanning from radio to 7— rays (9, 11) 
confirmed the position of the source to be coincidental with the nucleus of an inactive galaxy and 
showed the presence of rapid time variability 10 2 s) at high energies (8, 10) , indicative of a 
highly compact source of emission (~ 30 million km or ~ 0.2 AU). Furthermore, the detection 
of a relativistic and highly collimated radio jet, for the first time associated with a TDF (10, 11), 
makes Sw J1644+57 analogous to a small-scale blazar (10). Following the identification of 
this source as a tidal disruption candidate (TDC), we initiated a series of twelve, bi-weekly 
observations with the large X-ray satellite XMM-Newton, starting approximately 19 days after 
the BAT trigger (Fig. 1), together with a further Suzaku pointing ~ 10 days earlier (13). The 
long-term evolution in the X-ray (0.2-10keV; observed frame) luminosity followed roughly the 
expected rate of mass return to the black hole [oc t~ 5 / 3 for TDFs (2)] and the short-term, rapid 
variability was also apparent. 

We produced light curves for all observations of Sw J1644+57 over the 0.2-10.0 keV energy 
band and Fourier transformed these in order to obtain various power density spectra. A least- 
squares fits to the soft energy band (< 2keV) power spectrum for the Suzaku observation is 
consistent with a band-limited (red) noise component, described by a T —1.8 power-law 
plus Poisson (white) noise (Fig. SI and Table SI). However, the 2-10keV power spectra of 
both the Suzaku and first XMM-Newton observations (hereafter XMM #1) displayed a potential 
Quasi-Periodic Oscillation (QPO) component near 5mHz (Fig. 2). This feature has a centroid 
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frequency of v ~ 4.8 mHz and frequency width (full-width at half-maximum) of <5z/ suza k u < 
0.4 mHz and 5u xmm = 0.3 mHz (quality factor Q = u/5u > 12 and ~ 15 for Suzaku and 
XMM #1 respectively). The fractional root-mean- square (r.m.s.) variability in the QPO is > 
2.8% and ~ 4% for Suzaku and XMM #1, respectively (Fig. S2). 

Assuming the signal is on top of a purely Poisson-noise time series, the limit on the r.m.s. 
variability found here results (14) in a single-trial significance in the Gaussian limit of 3.8a for 
Suzaku and 2.2<r for XMM #1. However, in order to rigorously quantify the strength of the sig- 
nal outside of the Gaussian regime and fully account for the presence of missing and unequally 
spaced data (Fig. 1), we conducted Monte Carlo simulations. The method we adopted is based 
on well established procedures in timing studies of compact objects (14-16). A series of 50,000 
light-curves with the same average intensity, standard deviation, and number of bins as the orig- 
inal data were produced, and the noise power distribution was found at all Fourier frequencies 
and compared to the real observations (Fig. S4). In this manner, we found that the chance that 
the individual detections at ~ 5 mHz are due to random noise to be Pf a i se \ suzaku = 1-4 x 10~ 4 
and Pfaise\xmm < 5 x 10~ 4 for Suzaku and XMM #1 respectively (Fig. S7). After correcting for 
the initial "blind- search" of the Suzaku data and accounting for the fact that the QPO was de- 
tected in two independent observations, with different telescopes, the probability of two chance 
3a detections arising from random noise alone was found to be 1.52 x 10~ 5 . The observed quasi- 
periodic signal at ~ 5 mHz in Sw J 1644+57 is statistically highly significant (4.33<r assuming 
Gaussianity in the probability). 

QPOs are regularly seen in stellar mass black holes. Recently, QPOs have also been ob- 
served in a single AGN (17), and in a couple of potential intermediate-mass black holes (~ 
10 3 M ) (18, 19). Despite the lack of a unique physical explanation, most models (20-24) 
strongly link the origin of QPOs with orbits and/or resonances in the inner accretion disc close 
to the black hole. The detection of a QPO approximately 10 days after Sw J1644+57 became ac- 
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tive requires that an accretion disc was formed shortly after the start of the initial TDK The char- 
acteristic time (2, 25), £f a u, it takes for material from the disrupted star to fall back towards the 
black hole in Sw J1644+57 from a pericenter distance of R p ~ 13(M BH /10 6 M )" 5 / 6 R s (11) 
(where the Schwarzschild radius Rs = 2GM/ c 2 ) is tf a u 1 day. This is consistent with the for- 
mation of an accretion disc due to stream-stream collision (2) after a small multiple of the tf a n. 
The Keplerian frequency at the radius of the innermost stable circular orbit (ISCO = 0.5 — 3R S , 
depending on whether the black hole is spinning or not), just short of the black hole's event 
horizon, is generally the highest characteristic variability frequency. If the 5 mHz QPO centroid 
frequency is set by orbits at the ISCO, it would imply a black hole mass between ~ 4.5 x 10 5 
and 5 x 10 6 M - for a non-rotating and maximally rotating black hole, respectively. This is in 
line with predictions based on simultaneous X-ray and radio observations (26) (~ 3.2 x 10 5 M ) 
as well as the upper limit imposed by the M — L relation (< 2 x 1O 7 M ) (8, 10). Keplerian 
frequencies at the ISCO scale inversely with black hole mass; assuming a ~ 10 6 M black hole, 
the ~ 5 mHz QPO would correspond to a frequency of ~ 500 Hz around a ~ 10 M black hole. 
This is remarkably similar to the 450 Hz oscillation found for GRO J1655-40 (27). Our results 
thus confirms a fundamental aspect of disc/disruption theory (2) and highlights the scale invari- 
ant nature (28-31) of the underlying physics governing the accretion flow onto supermassive 
(~ 10 5-9 M ) and stellar-mass (~ 10 M ) black holes, several of which have displayed similar 
X-ray periodicities (32). 

Moreover, the Eddington limit for this black hole is < 6 x 10 44 erg s _1 , making the peak lu- 
minosity where the QPO is observed highly super-Eddington (Fig. 1). Thus, quasi-periodicities 
should be preserved in not only mildly super-Eddington accretion flows, as is the case for the 
only other supermassive black hole RE J 1034+396 (17,33) showing a QPO, and the stellar-mass 
black hole GRS 1915+105 (34), but also in potentially highly beamed, anisotropic sources. Re- 
cent numerical simulations have began to examine the role of relativistic jets in the production 
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of QPOs (24, 35), however, even in this scenario, it is the connection between the accretion disk 
and the base of the jet that give rise to the quasi- periodic signal. 

Although we have drawn from a rich literature that has ensured a high standard of robust 
techniques, there is still the concern that other uncertainties in the theoretical framework of ac- 
cretion flows could be at play and result in unaccounted systematic errors in the absolute noise 
model around black holes. Noise may not be purely white, purely red or a simple combination 
of the two, resulting in a smooth, featureless continuum. Indeed, (36), using general relativis- 
tic, magnetohydrodynamic simulations of the accretion flow onto a Schwarzschild black hole 
showed possible high frequency quasi-periodic features which were identified with properties 
of the turbulent flow. If this turns out to be the general case, QPOs similar to the one detected 
here would still be highly useful in determining the physics of the accretion flow, but might not 
be trivially related to the fundamental properties of the black hole, i.e. mass and spin. 
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Fig. 1. XMM-Newton and Suzaku light curve of Sw J 1644+57 together with the SWIFT-XKT 
0.3-10keV light curve for reference (grey) (37, 38). For twelve XMM-Newton observations 
(XMM#l- 1 2) we extracted 0.2-10 keV source and background light curves from the PN camera, 
using 40-arcsec circular regions. After accounting for the flaring background, a substantial 
fraction of the data in XMM #2, 3 and 12 had to be excluded. For Suzaku, we used a box 
region of 250-arcsec to extract the source light-curve from the two front-illuminated cameras 
and 150-arcsec for the background. Every observation had a background level significantly less 
than 5% of the total flux. For each observation we created an energy spectrum to which we 
fitted a model consisting of an absorbed power-law and used this to obtain a conversion factor 
between the count-rate and fluxes in physical units. The average flux levels in each observation 
are shown in real time in the main axis (stars, where the one s.d. error bars are smaller than the 
symbol size) and the inset compares each observation with data points binned in 32-s intervals. 
The vertical lines separate the various observations. The right-hand axis gives the conversion 
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Fig. 2a (Top): Power spectra for Suzaku and Fig. 2b (Bottom): for XMM #1, in the 2-10 keV 
energy range. The power spectra are normalised such that their integral gives the squared r.m.s. 
fractional variability. The Poisson noise level expected from the data errors is shown as the 
dashed horizontal line. We checked that the peak is robust to a variety of frequency and time 
resolutions. The arrow in both panels mark the presence of a QPO with a centroid frequency 
of Suzaku ~ 4.8 mHz. The solid curves enclose the range of the best fit to the underlying 
continuum as described in Table SI (Model 2). The dotted curves show the 99.99% and 99.73% 
(3a) threshold for significant detection. 
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1 Data Reduction 

1.1 Suzaku 

Suzaku observed Swift J164449.3+57345 1 (hereafter Sw J1644+57) on 201 1 April 6 for a total 
ontime exposure of 52 ks. The three operating detectors constituting the X-ray Imaging Spec- 
trometer (XIS (39)) on-board of Suzaku were operated in the normal, full-frame imaging mode 
with both front and back illuminated detectors in the 3x3 and 5x5 editing modes. Using the 
latest HEASOFT v6.1 1.1 software package we began by processing the unfiltered event files for 
each CCD following the Suzaku Data Reduction GuideQ Detailed good time intervals (GTIs) 
were produced using the FTOOL XISTIME and setting the option "bstgti=no". New attitude files 
were then created using the script AEATTCOR^] in order to correct for shift in the mean posi- 
tion of the source caused by the wobbling of the optical axis (40), however, we note that this 
effect was minimal and the correction was only used for consistency. The FTOOL XISCOORD 

1 http://heasarc.gsfc.nasa.gov/docs/suzaku/analysis/ 
2 http://space. mit.edu/cxc/software/suzaku/aeatt.html 

12 



was used to create new event files which were then further corrected by re-running the Suzaku 
pipeline with the latest calibration, as well as the associated screening criteria files. The good 
time intervals mentioned above were also employed in all cases to exclude any possible teleme- 
try saturations. XSELECT was used to extract spectral products from these event files with light 
curves being created with 16 seconds time resolutions - the minimum possible in this observing 
mode. The final good time exposure was 38.1 ks. 

In order to estimate the level of pile up suffered by the data we used the script (41) PILE J^ST^] 
to create a pileup map out of a Suzaku event data file. It was found that the data suffers very 
little from pile up (maximum pileup fraction of 5%) and therefore we chose to employ a square 
box region with a width of 240 pixels (~ 250") without excluding any data. Background 
spectra were extracted from a box region having a width of 150" elsewhere on the same chip. 
Individual ancillary response files (arfs) and redistribution matrix files (rmfs) were produced 
with the script XISRESfQ - which calls the tools XISRMFGEN and XISARFGEN - with the 
"medium" input. The resulting spectra had a 0.5-10 keV count rate of ~ 4.3 and ~ 4. lc/s for the 
front illuminated and back illuminated detector respectively, with the background constituting 
~ 0.4% and ~ 0.6% of the total flux. 

1.2 XMM-Newton 

XMM-Newton performed twelve, approximately bi-weekly pointing observations of Sw 11644+57 
starting on 201 1 April 16 and finishing 201 1 October 3, with each observation lasting between 
~ 16 ks to ~ 26 ks. In all cases, the EPIC-pn camera (42) on-board XMM-Newton was operated 
in "imaging" mode with a "medium" optical blocking filter. The EPIC-MOSl and EPIC-MOS2 
cameras (43) were operated in the "imaging" data mode with the EPIC-MOSl and EPIC-MOS2 
being in the "PrimePartialW2" and "FastUncompressed" submodes respectively. In this work 

3 http://space. mit.edu/cxc/software/suzaku/pest. htm 
4 http://suzaku.gsfc. nasa.gov/docs/suzaku/analysis/xisresp 
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we only use the EPIC-pn data, from which we always start with the unscreened level 1 data 
files. We generate concatenated and calibrated event lists using the latest XMM-Newton Science 
Analysis System v 11.0.0 (SAS). EPIC-pn events were extracted from a circular region with ap- 
proximately 40" but varying slightly depending on how far the observation is from the chip gap. 
Bad pixels and events too close to chip edges were ignored by requiring "FLAG = 0" and "PAT- 
TERN < 4". Flaring background meant that a significant fraction of of the data in XMM #2, 
3, and 12 had to be discarded. The energy channels in all observations were initially binned 
by a factor of five to create an energy spectrum and both source background light curves were 
produced in the 0.2-10keV energy range. We used the SAS task EPATPLOT to assess the level 
of pile-up in the EPIC-pn spectrum and found it to be insignificant. The highest EPIC-pn count 
rate is approximately 10 count s -1 for XMM #1, which is below the nominal pile-up limit for 
the instrument in "LargeWindow" submode. Response files were created in the standard way 
using the tools RMFGEN and ARFGEN. The total good exposure time for XMM #1 was 16.96 ks. 
Finally we rebinned the energy spectrum with the tool PHARBN^J to have 3 energy channels per 
resolution element, and at least 20 counts per energy channel. 

2 Energy Spectra 

2.1 Suzaku 

We analysed briefly the energy spectra within the X-ray fitting package XSPEC (44). A model 
consisting of an absorbed power law (ZPHABS xZPOWERLAW in XSPEC ) resulted in a fit with 
X 2 /D.o.F = 1181.2/779 (Reduced-^ 2 = 1.516) which is not formally acceptable. However, as 
can be seen from Fig. Sla (top and middle panels), quite a number of the systematic residuals 
are in fact coming from the region between ~ 1.6 — 2.0 keV, likely due to the known systematics 
uncertainties in the effective area near the silicon K and gold-M edge. Excluding this region 
5 http://virgo.bitp.kiev.ua/docs/xmm_sas/Pawel/reduction/pharbn 

14 



Suzaku XMM #1 
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Fig. Sla (left): Front (black) and back (red) illuminated Suzaku spectra (top) fit with an ab- 
sorbed powerlaw. The middle panel shows the residual to this fit over the full energy range and 
the bottom panel after excluding the silicon region, the front illuminated data below 1 keV and 
adding a phenomenological disk black body. Fig. Sib ( right) Same for the XMM #1 EPIC-pn 
spectrum in the full energy range. Note the severe drop in total counts below « 1 keV as well 
as the clear presence of un-calibrated silicon features just below 2 keV in the Suzaku data. 

improves the fit to x 2 /D.o.F = 982.6/723 (Reduced-^ 2 = 1.36) with now the only noticeable 
systematic residues coming from the front illuminated data below ~ 1 keV as a result of its 
smaller effective area as compared to the back illuminated camera. Excluding the FI data below 
1 keV results in x 2 /D.o.F = 833.9/688 (Reduced-^ 2 = 1.21). 

SWIFT observations of Sw J1644+57 during the first two days after its detection showed that 
the system exhibits times where the spectrum is well fit by an absorbed powerlaw at comparably 
low-flux states and preferably a black body at higher-flux states (45). In fact, the authors suggest 
that during those observations, the spectra are better fit by a superposition of a thermal disk with 
a temperature of ~ 1 keV, and a non-thermal component. A similar study using the combined 
SWIFT-XKT data for the first 10 days after the trigger (8) finds that the spectra at various 
intensity levels are comparably well fit by either an absorbed log-parabola model, absorbed 
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broken-powerlaw model or an absorbed powerlaw model plus an extra soft thermal component. 
However, a simple absorbed powerlaw model was inferior to the other combination during the 
first 10 days. The authors also find that during that time, the spectra hardens from T ~ 2.3 
to T ~ 1.8 as the intensity increases. A further study (9) also explore the SWIFT-XKT data 
during the first ~ 15 days and find a time-averaged photon index of T = 1.80 ± 0.25 when 
modelling with an simple absorbed powerlaw, but also states that a powerlaw with an index of 
r = 1.6 having an exponential cutoff at ~1 keV can also model the data equally well. This latter 
phenomenological combination has the same effect as an extra soft component affecting the 
spectra below ~ 1 keV, similar to the phenomenological addition of a thermal disk component 
used by (8,45). 

We explore whether the addition of a DISKBB component significantly improve the quality 
of the fit and find that it provides only a mild improvement of A% 2 = 27.7 for 2 degrees of free- 
dom. This fit (shown in the bottom panel of Fig. Sla) has an absorbing column iV H ol^](1.81 ± 
0.06) x 10 22 cm~ 2 , a powerlaw index of 1.83 ± 0.02 and a disk temperature of 0.35lo;og keV. 
The unabsorbed 0.5-10keV powerlaw and DISKBB fluxes are 1.93lo;° 2 x 10~ 10 erg cm -2 s _1 
and 8^ x 10~ 12 erg cm -2 s -1 , respectively. 

It should be stressed that it is often the case with active galactic nuclei that energies below a 
few keV show the presence of a broad and smooth "soft excess" (see e.g. (46)) which, despite 
the observed AGN hosting black holes spanning many orders of magnitude in mass, remains 
rather constant in temperature when fitted with a phenomenological disk black body (e.g. (47)) 
suggesting that the excess is in fact atomic in nature and not the true disk emission which often 
peaks in the UV. As such, the temperature quoted above is only a fiducial disk temperature and 
should not be though of as the actual accretion disk temperature. 

6 A11 errors shown in this section are 90 per cent confidence errors obtained by allowing all parameters to vary. 
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2.2 XMM-Newton 

We find that a model consisting of an absorbed powerlaw (ZPHABS x ZPOWERLAW in XSPEC ) 
provides a formally acceptable fit to all twelve XMM-Newton spectra, with reduced-^ 2 ranging 
from 0.9 to 1.14. Simultaneous fit to all 12 spectra, allowing for variation in all model pa- 
rameters between observations, resulted in x 2 /D.o.F = 8190.3/8219 (Reduced-^ 2 = 0.997). 
The data shows a distinct hardening of the powerlaw index with time with the power law index 
going from 1.72 ± 0.01 in the first observation (day ~ 19) to 1.52 ± 0.02 for XMM #2 (day 
~ 33 and settling at approximately 1.3 at later observations (for XMM #6-12, T < 1.4). Figure 
1 in the main journal shows the evolution in the flux level of the twelve XMM-Newton obser- 
vations, together with the preceding Suzaku pointing as well as the SWIFT-XKY light curve for 
reference. 

The total column density varies in the range ~ (1.32 — 1.62) x 10 22 cm~ 2 and assuming 
the spectra are indeed simple powerlaws, we find that this variation is statistically significant, 
in contrast with (9). A model having a static value for Nu of (1.41 ± 0.01) x 10 22 cm" 2 yields 
X 2 /D.o.F = 8441.9/8230 (Reduced-^ 2 = 1.03). However, if we add a phenomenological 
soft-excess in the form of a DISKBB, an equally valid fit (x 2 /D.o.F = 8175.8/8206; (Reduced- 
X 2 = 0.996) can be obtained with a static Nn and variable disk flux and temperature. In 
the first five XMM-Newton observations where the temperature of this fiducial soft-excess is 
well constrained^} we find its value to range between 0.24-0.36 keV (compare with the Suzaku 
value of 0.35+o;ogkeV). Fig. Sib shows the EPIC-yn spectrum of XMM #1 fit with a single 
powerlaw (top and middle panels) as well as after the inclusion of a DISKBB (bottom panel). 
The smoother decline in the soft flux of XMM #1 as compared to Suzaku is likely due to the 
lower Nn of (1.36 ± 0.04) x 10 22 cm" 2 compared to (1.81 ± 0.06) x 10 22 cm" 2 for Suzaku. 
7 Defined here as having errors less than 50% its value. 
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Fig. S2a (left): Average power spectrum for Suzaku and for XMM #1 Fig. S2b (right) obtained 
for energies below 2 keV. Note the differing frequency ranges. 

3 Power Spectrum 

3.1 Suzaku 

Suzaku's low orbit means that most targets are occulted by the earth for about a third of its 
~ 5.8 ks orbit. For this reason, the 38.1 ks observation in fact consists of fourteen periods 
ranging from ~ 1.98 ks to ~ 3.0 ks. To create a combined power spectrum, we broke the full 
observation into these fourteen continuous time series each having identical frequency bins and 
averaged their power spectra. We used 256 s intervals to produce a total of 99 frequency bins 
ranging from 2.44 x 10~ 4 Hz to 3.125 x 1CT 2 Hz with frequency resolution of 1.22 x 10~ 4 Hz 
up to 1.68 x 1CT 2 Hz and twice that at higher frequencies (geometrically binned) in order to 
produce the individual power spectra. A visual inspection of the averaged 0.2-2.0 keV power 
spectrum displayed no distinct, narrow features with a model consisting of Poisson noise and 
a power law component yielding a best-fit (least-squares method) of \ 2 l v — 115.2/96 and 
r = -1.82 ± 0.19 (Ax 2 = 2.71 for 1 D.o.F). This is shown in Fig. S2 (left). However, 
the spectrum above 2 keV (Figure 2 of the main journal) shows a significant feature at around 
5mHz. Visual inspections of power spectra having twice and four times the frequency resolution 
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Table SI: Fits to the power spectra of Suzaku and XMM #1. 

Model 1 Model 2 Model 3 





Suzaku 


XMM#l 


Suzaku 


XMM#1 


Suzaku 


XMM#l 


Poisson level 


0.30 ±0.02 


0.872 ±0.014 


0.32 ±0.02 


0.872 ±0.014 


n 09+0.02 

u ' oz -0.01 


0.871 ±0.014 




-1.79 ±0.12 


1+0.4 

^■■■-o.e 


i q+0.4 

^-o.s 


-1(f) 


1 q+0.4 

i -°-0.3 


-1(f) 


r 2 






-2.90±g;g2 


-3(f) 


q 1O + 0.55 

°- lo -0.03 


-3(f) 


Fbreak (mHz) 






i-4±S:S 


< 1.7 


1 4+0-4 


< 1.7 


iVi (xicr 5 ) 




0.095 ±0.037 


75 ± 15 


330 ± 140 


75 ± 15 


330 ± 140 


iV 2 (xlO- 9 ) 






sits 


7±4 




6±4 


v (mHz) 










a 01 +0.08 
4 ' 8i -0.03 


4 7+0-5 
+•'-0.2 


i5^(mHz) 










< 0.4 


298+ ' 001 


N v 










3.1 ±1.5 


3.6 ±2.0 


X 2 /" 


156.1/96 


392.5/250 


136.0/93 


384.9/249 


125.6/90 


375.9/246 



Note - Model 1 consists of a power law together with Poisson noise. Model 2 introduces a break to the power law 
and Model 3 adds a Lorentzian component to account for the presence of the QPO. In the case of XMM #1 , the two 
power law indices could not be constrained so we froze them at the values shown. 

also indicates the presence of this feature. A fit with the same model as to that for the spectrum 
below 2keV, i.e. a constant together with a power law, results in a fit with x 2 / u — 156.1/96 
(Model 1; Table SI) and an index of V = —1.79 ±0.12. Adding a break to this power law so 
that below a frequency F &reafc the continuum is described as Ti and above as T 2 resulted in a 
significant improvement to the continuum ( (x 2 /v = 136.0/93; Model 2; Table SI). From Fig. 
S6 we can see the visual improvement of a broken power law over the simpler continuum. 

An equally satisfactory description of the continuum (x 2 /v = 140.6.1/94) can be achieved 
by replacing the broken power law with the sum of two broad Lorentzian components (48-50), 
one having a centroid frequency frozen at zero and the other at a frequency similar to the break 
frequency above. In all cases, however, the presence of a feature above the continuum around 
5mHz remains and we finally model this QPO with a Lorentzian having a centroid frequency of 
4.8llg;o3 mHz and a quality factor Qsuzaku = Su/u > 12. 



19 



3.2 XMM-Newton 

We produced power spectra for all XMM-Newton observations, and the resultant power spec- 
tra were averaged together. The presence of a feature at a similar frequency to that shown for 
Suzaku above was immediately clear; however, by looking at the individual power spectra from 
each observation we found that the signal was driven by the first pointing. This is not surprising 
given the clear change in the spectral properties of the source between XMM #1 and the latter 
observations mentioned above, and the closer resemblance XMM #1 has with the spectral shape 
of (and proximity in time with) the Suzaku observation. Indeed, QPOs are known to be transient 
features with a good example being the QPO found in the active galaxy RE J 1034+396 where 
the 5.6(7 discovery (17) was not seen in three follow up observations nor in a preceding ob- 
servation (51). As such, we concentrate below on the description of the power spectrum made 
for XMM #1. Figure S2b shows a power spectrum for energies below 2keV fit with a power 
law with index T = —2. 261^2 (X 2 /^ — 341.6/250). In contrast to Suzaku, the corresponding 
power spectrum shows visual sign of a feature at ~ 5mHz even at this soft energy range. 

The lack of any significant feature below 2 keV in the Suzaku data can be understood in the 
context of the previous section in that not only does the total counts for Suzaku drops dramat- 
ically at these energies, it also suffers from clear calibration problems and confusion between 
possible variation in the column density or the presence of a further soft-excess component. As 
is clear from Fig. SI, the spectrum of XMM #1 does not contain a comparatively strong drop 
in counts, partly due to its higher effective area as well as possibly lower column density or 
stronger soft-excess component at low energies, nor does the spectrum suffer from distinct cal- 
ibration. These various factors propagates into a detection of a feature in the power spectrum 
below2keV in XMM #1. 

As detailed in the previous section, we cannot break the degeneracy found at the low en- 
ergy range between a model having a variable column density with a powerlaw continuum or a 
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static column with a further variable soft-excess. With this in mind, as well as to be consistent 
with Suzaku where a further difficulty arises due to low effective area and calibration uncertain- 
ties, we proceed through the remainder of the manuscript by examining the data above 2 keV. 
Figure 2 in the main journal shows the power spectrum at harder energies similar to that of 
Suzaku. We used 2048 s intervals to produce a total of eleven power spectra which were aver- 
aged to make a single power spectrum with 253 frequency bins ranging from 4.88 x 1CT 4 Hz to 
4.998 x 10 _1 Hz with frequency resolution of 2.44 x 1(T 4 Hz up to 3.37 x 1(T 2 Hz and increas- 
ing at higher frequencies (geometrically binned) in order to produce the individual power spec- 
trum. A visual inspection of the averaged power spectrum shown in Figure 2 of the main journal 
shows a significant feature at around 5mHz. Fitting a single power law to the continuum yields 
r = -2.2ll°; 4 with X 2 /v = 392.5/250. Replacing the single power law with a model similar 
to that of Suzaku, having Ti and T 2 frozen at -1 and -3 respectively results in an improvement to 
the continuum with A% 2 /Az/ = 7.6/1 (Table SI). A similar best-fit is obtained when replacing 
the broken power law with a combination of two broad Lorentzian (x 2 /v = 385.4/248), again 
in line with the results presented for Suzaku above. Adding a Lorentzian to model the QPO at 
~ 5mHz, we find that the QPO in the XMM #1 data has a centroid frequency of 4.7+Q 2 mHz 
and a quality factor Qxmm ~ 15. Within the precision of the measurements, the frequency of 
the QPO appears to be stable between both independent observations. 

3.3 SWIFT-XRT Timing Analysis 

Early results (8) claimed to find no evidence for any statistically significant QPO in the SWIFT- 
XRT curve during the first twenty days of the outburst. We have performed an analogous Fourier 
analyses on the data belonging to the first 20 days of the observations and summarise the results 
in Fig. S3. Due to the random observing cadence employed by the XRT, the power spectrum 
is highly aliased and dominated by noise making the various apparent features not statistically 
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Fig. S3a (left): Swift- XRT power spectrum made from the combined data from the first ~ 20 
days after the outburst. The presence of various features at low frequencies are immediately 
obvious and are likely a result of aliasing due to the random observing cadence employed by 
the XRT. Fig. S3b (right) Close up at higher frequencies showing the possible presence of 
a feature clearly above the local noise level, which, when fit with a Lorentzian, results in a 
centroid frequency of ~ 6.2mHz and a quality factor Q « 20. We emphasised that due to 
the noise characteristics of this power spectrum, this feature is not statistically significant. The 
frequency of the QPO reported in this manuscript is marked by the vertical arrow. 

significant, in agreement with early results. Interestingly, although not statistically significant, 
if we zoom into a narrower frequency range centred at ~ 5mHz (Fig. S3b) we do see a potential 
feature which, when modelled with a Lorentzian on top of Poisson and red noise, results in a 
centroid frequency of ~ 6.2mHz and a quality factor Q 20. 

4 Monte Carlo Sampling 

In order to confirm the accuracy of the estimated values presented in Table S 1 , we proceeded by 
generating a sample of 10,000 light curves sampled by the empirical distribution of the observed 
light curves themselves. Let q be the count rate and c ee r ; i the associated error of the i th time 
bin in the observed light curve. Based on this input light curve, 10,000 bootstrapped (52) light 
curves were simulated using the following prescription: the i th time-bin of a given simulated 
light curve have a count-rate of ccj and count-rate error c err;i with cq being drawn from a 

22 



normal distribution N(ci, c crr -i). The GNU Scientific Library (GSL) function gsLran_gaussian 
was used to draw a Gaussian random variate with mean zero and standard deviation c eir;i . The 
value returned by gsLran.gaussian was added to q to obtain cq. Since observed fluxes cannot 
be negative, we set cq to be equal to zero when its value returned by the above prescription was 
<= 0. A similar approach was used to generate GRB light curves by Stamatikos et al (53) (see 
also (54)). 

For each light curve, we generated power spectra in an identical manner to that described 
in the previous sections and modelled the data with the template described in Model 3, with 
all parameters except for the centroid frequency of the Lorentzian describing the QPO free to 
vary. Figure 2 shows the fractional root-mean-square (r.m.s) variability in the QPO as found by 
the integral under the Lorentzian in the "r.m.s. normalised" power spectra, ie. it is the square 
root of the product ttN^Su/2. The QPO in the Suzaku data is clearly very narrow (See Figure 2 
in the main journal), possibly lying within just one frequency bin. The upper limit shown in 
Table SI give rise to the low peak in Fig. S4 at ps 2.8%. We checked this limit on the fractional 
variability by further running a series of simulations with the width of the Lorentzian frozen at 
4 x 1CT 4 , 1 x 1CT 5 and 5 x 1CT 6 Hz, and in all cases we obtained a lower limit consistent with 
the result shown in this figure. For XMM#l, we find that the r.m.s. variability in the Lorentzian 
is better constrained at p=* 4% compared to that of Suzaku, in line with the constraints on the 
width and normalization quoted in Table SI. 

Assuming that the signal is on top of a purely Poisson-noise time series of length T and 
count-rate I x , which is the product of M individual power spectra, the limit on the r.m.s. vari- 
ability (r) found here, combined with Equation 2.17 of (14), results in a single-trial significance 
in the Gaussian limit > 3.8a for Suzaku and > 2.2a for XMM #1. However, in order to rigor- 
ously quantify the strength of this QPO outside of the Gaussian and Poissonian-noise limit as 
well as to account for all the frequencies probed in the various power spectra, we perform in the 
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Fig. S4a (left): Fractional root-mean-square variability in the Lorentzian component modelling 
the QPO in 10,000 bootstrap samples for Suzaku. We generated a large number of light curves 
by sampling the mean and standard deviation in each time step of the original light curves for 
both Suzaku and XMM#\ (Fig. S4b (right)). For each of the new light curves we then generated 
a power spectrum in an identical manner to the observed data and fit model 3 (Table SI). From 
the normalisation and equivalent width of the Lorentzian component used to model the QPO, it 
is then trivial to obtain the fractional r.m.s. variability and an estimate of the error in this value. 

following section a full series of Monte Carlo simulations. 

Fig. S5 shows the distribution in the poweif^present in the QPO bin with centroid frequency 
of 4.88 mHz after having the continuum-model level subtracted (i.e. (I — M) f =4: 88mHz ). The 
continuum in each of the 10,000 simulated power spectra was found by fitting Model 2 (Ta- 
ble SI) to each power spectrum. The dashed red line shows that the observed value lies at the 
peak of the distribution in both observations, which confirms that the signal level in the QPO 
bin are not a statistical fluctuation in either observations. 

5 Monte Carlo Simulations 



The previous tests and quoted significances accounts for white noise variability; however, the 
presence of red noise is known to be a potential source of false features in the X-ray power spec- 



8 Throughout the rest of this work, use power, 7/ to refer to the signal strength in any frequency bin, / when 
the power spectra is normalised such that their integral gives the squared r.m.s. fractional variability. 
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Fig. S5a (left): Histogram of the power in the 4.88 mHz frequency bin in 10,000 bootstrap 
simulations for Suzaku. The power (minus continuum level) in the frequency bin with centroid 
frequency of 4.88 mHz in the 10,000 sampled spectra are shown as a blue histogram and it 
can be seen that it agrees well with the real observed value (red vertical dashed lines). This 
addresses the robustness of the signal in the QPO bin against possible statistical fluctuations. 
Fig. S5b (right): shows the same forXMM#l. 

tra of X-ray binaries (55). In order to fully investigate the significance of the aforementioned 
feature, we implemented Monte Carlo techniques and compared the spread of 50,000 simulated 
power spectra to the observed power spectra in order to estimate the significance, in a similar 
manner to that outlined in Uttley, McHardy & Papadakis (16). 

We simulated 50,000 light curves using the method of Timmer & Koenig (5(5), for a broken 
power law continuum similar to that of Model 3 (Table SI), having the low frequency index 
and break frequency frozen at their respective value. The high frequency power law slope was 
incremented between -2.4 and -3.2 in steps of 0.1 and the simulated light curve was resampled 
and binned to match exactly the pattern of the observed light curve. Power spectra for all 
simulated light curves were produced in an identical manner to that described in the previous 
sections. 

Fig. S6 (top) shows a characteristic simulated power spectrum (red) having continuum 
parameters similar to that of Model 3. Also shown is the real data (black) where we can see that 



25 



1000 



o 

Q. 

o 

N 



o 
Z 



100 



10 



1000 



& 
o 
a. 

o 

N 



o 
Z 



100 



10 



Suzaku 



Suzaku observation 
Sample MC bkn-power 




Suzaku observation 
Sample MC unbroken 




100 



o 
a. 

o 

N 



o 
Z 



10 



100 



o 
-■ 

N 



o 
Z 



10 



XMM#1 



XMM #1 observation 
Sample MC bkn-power 




Poisson noise 




XMM #1 observation 
Sample MC unbroken 




10~ 



0.01 

Frequency (Hz) 



0.1 



Frequency (Hz) 

Fig. S6 Sample simulated power spectrum from a broken (top) and unbroken (bottom) power 
law continuum. The broken power law provides a better description of both Suzaku (left) and 
XMM #1 (right) data (shown in black; see also Table SI). One of our statistical tests involves 
simulating 50,000 such spectra and comparing the results to the real data (see Figs. S7-S9). 



the simulated spectra indeed display similar noise level and characteristics as the original data. 
The bottom panels display a similar spectrum simulated for an unbroken power law (blue). For 
each of the 50,000 Monte Carlo power spectrum simulated as described above, we used Model 2 
(Table SI) to find the best fit to the red plus Poisson noise continuum. We show in Fig. S7 the 
data-to-model ratio for the observed data (solid histograms) above the low frequency break as 
well as for the mean of all the simulations (blue dotted lines). We used Monte Carlo simulations 
to estimate 1 s.d errors in each frequency channel as exemplified in Fig. S8 and show in Fig. 
S7 the 3cr (99.73%) and 99.99% confidence limits for the power as a function of frequency in 
the Monte Carlo simulation (dashed curves) . 

As exemplified in Fig. S8, for each simulated power spectrum, we found the spread in the 
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Fig. S7a (left): Data- to-model ratio for Suzaku. Fig. S7b (right) similarly for XMM #1. The 
ratio between the observed power spectra and the continuum model (Model 2) are shown as the 
solid histogram. The dotted curve shows a similar ratio for the average of 50,000 Monte Carlo 
realisations to their respective continuum model and the dashed curves shows the 3cr (99.73%) 
and 99.99% confidence limits derived from the standard error obtained from the 50,000 realisa- 
tions in each frequency channel. The arrow in both figures indicates the centroid frequency of 
the QPO as found in the Suzaku data. It is clear that the XMM #1 observation closely matches 
the frequency in the Suzaku. 

value (I/M)f and obtained the spurious-detection-probability, -P/ a ; se , that this value would be 
at or greater than the value found in the real data at the frequency bin of the QPO. Fig. S8 
shows the distribution for (I/M) f=4.88mHz with the observed value highlighted in red, and Fig. 
S9 summarises this result for all Fourier frequencies above the low frequency break. For 50,000 
simulations, we found a power above the continuum at ~ 4.88 mHz similar to that of the real 
data only 7 times for Suzaku. As hinted in Fig. S7, the QPO in the XMM #1 observation, despite 
being at the same frequency as that found for Suzaku, in fact occupies multiple frequency bins. 
This can be easily seen in Fig. S7, where there are two neighbouring bins that reaches the 
99.73% curve. We show in Fig. S9 that the chances of any two (or 3) neighbouring bins having 
a similar integrated power as the observed QPO bin is indeed very low at < 5 x 10~ 4 . 

So far we have calculated the probability of finding a QPO with the same power level as 
that found in the real data in a given frequency bin (i.e. Figures S8 and S9). Since in general, 
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when we first inspected the Suzaku power spectrum, a QPO would have been reported at any 
plausible frequency we proceed by calculating the probability of chance occurrence e = 1 — 
rin=X.in Pfaise\n), where P false is the false detection probability from Fig. S9, of a detection 
at any Fourier frequency greater than f min which we have set to be the low frequency break for 
convenience. This product effectively accounts for our initial "blind- search" of the Suzaku data 
and reduces the single observation significance to e = 2.31 x 1CT 2 (2.27cr). By itself, the Suzaku 
observation would not constitute a discovery. 

All our work so far however, has only considered significance in individual observations. In 
both cases that we see the QPO, our extensive set of Monte Carlo simulations has shown that for 
any reasonable choice of binning and assumptions the single trial significance at the frequency 
of the QPO (~ 5 mHz) is greater than 3a. This is highlighted in Fig. S9. We have shown 
here that the QPO exits in data based on two different telescopes which of course significantly 
reduces the likelihood that it is due to instrumental systematic effects. We have seen it in the 
two spectra with similar flux levels and spectral shape and at the same frequency. The latter is 
particularly important because when we embarked upon "searching" for a QPO in the XMM #1 
data, we did not do so "blindly". Combining the blind-chance probability of a QPO in the 
Suzaku data with the two-trials^] probability for a QPO at the same frequency in the XMM #1 
data, 6.6 x 10~ 4 (3.41cr), we find the probability that these two detections arises from random 
noise alone is 1.52 x 1CT 5 , which assuming a Gaussianity equates to a detection at the 4.33a 
leve{j^] Thus the many statistical tests and simulations makes it clear that the QPO is robust and 
statistically significant. 



9 Using the 2 bins spurious detection probability in Fig. S9 and correcting for the fact that we effectively 
searched two frequency bins and not one. 

10 Had we instead used the 3 -neighbouring bins probability, accounting for the extra number of trial, we would 
have detection at the 4.43c level of confidence. 
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Fig. S8: Histograms (top) and cumulative frequency (bottom) for the power in the QPO bin. The 
ratio of the power in the frequency bin with centroid frequency of 4.88 mHz to the continuum 
model in the various Monte Carlo simulated spectra are shown as a histogram and cumulative 
frequency for Suzaku (left) and XMM #1 (right). For comparison, we show the real observed 
value as the red vertical dashed lines. The chance of making a spurious detection (the false alarm 
probability) at the same power level as the Suzaku (XMM #1) observed data at this frequency is 
1.4 x 10~ 4 (6.8 x 1CT 3 ) based on a total of 50,000 Monte Carlo simulated power spectra similar 
to that shown in Fig. S6 (top). Note however that this is the probability for a single bin which, 
although appropriate for Suzaku, is not the case for XMM #1 where the QPO is clearly seen in 
at least two frequency bins (see Figs. S7 and S9). 
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Fig. S9a (left): Spurious detection probability for Suzaku. For all frequency channels higher 
than the low frequency break (Fig. S6), we found the probability that the ratio between the 
power and continuum at the specific frequency probed would be greater than the level found for 
the QPO in the real observation. The vertical dashed line shows the centroid of the frequency bin 
containing the QPO and the horizontal dotted line shows the sensitive limit due to our limited 
number of simulations (50,000). The QPO is within a single bin in the Suzaku observation (see 
Fig. S7) and we recover the value shown in Fig. S8 for Suzaku at ~ 5 mHz (1.4 x 10~ 4 ). 
Fig. S9b (right): Shows a similar analyses for XMM #1. However, it is clear from Fig. S7 that 
the QPO in the XMM #1 observation consists of multiple bins (at least 2) and indeed integrating 
over neighbouring bins we find that the chance of any two or three neighbouring bins having the 
same integrated power as the QPO is significantly low at < 5 x 10~ 4 . The significance levels 
are calculated assuming a Gaussian distribution from the percentage level shown on the left of 
the horizontal dashed lines. 
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